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1. Daubechies wavelets 

The polynomial formed by the first N terms of the binomial expansion 

k=0 ^ " ' k=0 

that is, the polynomial 

n , , ^ fk + N -1\ k at N(N + 1) 2 {2N-2\ jv— i 

p N (y) = J2( k )y k = i + Ny + ^- 2 — -y +■■■ + [ N _ 1 jy . a- 1 ) 

k=o v ' v ' 

where (a)f. = T(a+k)/T(a), plays an important role in the construction of the compactly 
supported Daubechies wavelets. There is a close connection between the zeros of Pn(v) 
and the 2iV filter coefficients h(n) of the Daubechies wavelet D 2 n- For a complete 
account of the theory we refer to Chapter 6 of Daubechies (1992). In this section 
we give the details that are relevant to this paper. In later sections we describe the 
asymptotic methods for obtaining the zeros of Pn(v) and discuss methods how to obtain 
the coefficients h(n) numerically. This paper is an extension of work done by Shen &: 
Strang (1995). 
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1.1. Some properties of Daubechies wavelets 

We recall that the filter coefficients h(n) define the dilation equation 



2N-1 



(f)(x) = V2 h(n)^>(2x ~ n ), (1-2) 



n = 



the solution of which is called the scaling function. We take the following normalization 
of 4> and, hence, of the coefficients h(n): 

oo 2N-1 

/ <f>(x)dx = l, J2h(n) = V2. (1.3) 

For the Daubechies wavelets the filter coefficients are real. 

When the filter coefficients and the scaling function <j> are available, the correspond- 
ing compactly supported orthogonal Daubechies wavelet is given by 

i 

ij){x) = y/2 Y (-l) n h(l-n)</>(2x-n), (1.4) 

n=2-2N 

which is also denoted by D 2 n- 

Denote the Fourier transform of <f> by 

<(>{£) = -— l e- ix ^(x)dx. (1.5) 
V 2ir J-oo 

Then (1.2) and (1.4) give the relations 

HO = m {i/2) £(£/2), fa) = M o(C/2) 0(€/2), (1.6) 



where 



2iV-l 

™o(0 = 7f E Hn)e~^, 



n=0 
1 



/*o(0 = -^ E (-l) B Ml-n)c- faC 

^ 2 n=2-2W 

= -e _i€ m (-£-7r). 



(1.7) 



On the one hand, the orthonormality of the functions <j>{- + k),k € 7L can be 
expressed in terms of the filter coefficients by the relation 

J2Hn-2j)h(n-2k) = S jik , 

where j,k,n run over all relevant indices that define the finite set h(n). On the other 
hand, the orthonormality is given by the relation 

ho(C)| 2 + |m (C + 7r)| 2 = l. (1.8) 
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The trigonometric polynomial mo(£) plays an important role in wavelet theory. It 
is called the filter function or transfer function. In classifying the compactly supported 
orthogonal wavelets, Daubechies assumes that the mo(£) has an N— fold zero at ±7r. 
That is, she writes 

mo(0=(—^—) ^/(n)e-^. (1.9) 

A consequence is the N— fold zero of m (£ + -k) and m (— £ — -k) at £ = 0. Using the 
relation for ip in (1.6) and fio in (1.7) we see that the N— fold zero of mo(£) at — n is 
related to the vanishing of the N moments 

J ip(x)x k dx, k = 0, 1 . . . , N - 1. 

This property gives a certain degree of approximating quality of the wavelet and leads 
to vanishing moments for the filter coefficients: 

2JV-1 

(-l) n n k h(n) = 0, k = 0,1,... N- 1. (1.10) 

n = 

1.2. The construction of Daubechies wavelets 

Taking the symmetric product of the two filter functions, and writing z = e^, we have 
m (0m (-0 = \mo(0\ 2 = cos 2W (^) i^fo), (1.11) 
where Pn(v) is the polynomial defined by 

N-l 

Pn{v) = Qn(z)Qn(1/z), Q n (z) = /H^" n (1-12) 

n = 

and y is defined by z + 1/z = 2 — 4y, that is, y = sin 2 Indeed, the product 
Qn(z)Qn(1/z) can be written as a polynomial of y or cos£, as can be verified easily. 

When the coefficients f(n) in (1.12) are known, the filter coefficients h(n) follow 
from (1.9) and (1.7). The problem for constructing Daubechies wavelets is to find f(n). 

Substituting (1.11) into the orthonormality relation (1.8) we obtain for the poly- 
nomial Pn the relation 

(l-y)"P N (y) + y"P N (l-y) = l. (1.13) 

A polynomial solution of this equation is the polynomial introduced in (1.1). To see 
this we introduce the incomplete beta function 

7 *M) = W~JS f ^{l-tf-Ut, B(a,b) = ^^-, (1.14) 
B(a,b) J T(a + b) 
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where > 0, 5R& > and B(a,b) is the (complete) beta integral. We have I x (a,b) + 
Ii- X (b, a) = 1 and we verify easily that the function 

P N (y) = (l-y)- N I 1 . y (N,N) (1.15) 

satisfies the relation in (1.13). Also, this Pn is the polynomial introduced in (1.1). 
Namely, using the substitution t = (1 — y)u we obtain from (1.14): 

^/iV-l\ (iV + fc-l)!(iV-^-l)! 

h \ k J (N-iy.(N-iy. 




1.3. Determining the filter coefficients 

When we know the quantity P^r, we can try to solve for Qat hi (1.12), and/or to find 
the f(n). As an example consider N = 2. From (1.12) we obtain 

l + 2y = [f(0) + f(l)/z}[f(0) + f(l)z] 

= /(0) 2 + /(l) 2 + /(0)/(l)(z + l/z) 
= ([/(0) + /(l)] 2 -4 2/ /(0)/(l), 

from which we can solve for /(0),/(l). From (1.3), (1.7) and (1.9) we already know 
that /(0) + f(l) = 1. Furthermore we obtain [/(0) - /(l)] 2 = 3. This gives /(0) = 
(l + V / 3)/2,/(l) = (l--\/3)/2. The filter coefficients follow from the relation (see ((1.7) 
and (1.9)) 

[(1 + l/,)/2] 2 [/(0) + /(l)/z] = -L[h(0) + + M2)A 2 + H3)/z% 

This gives the coefficients of D 4 : 

1 + ^ 3 + ^3 3- V3 1- V3 

M0) = ^7T' M1) = ^7T' M2) = ^7T' M3) = ^7T- (L16) 

For iV = 3 we can also obtain exact solutions of the nonlinear equations, but for larger 
values of N this is not possible. In addition, the complexity of the computational scheme 
increases. 
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Remark 1.1. The above equations for /(0), /(l) are symmetric with respect to these 
quantities. Interchanging the values of /(0),/(l) gives a different set h(n) (with h(j) — > 
M3 = 0,1,2,3). The present choice of /(0),/(l) gives Q 2 (z) = f(0) + f(l)/z, 

with a zero at z\ = 2 — y/3 , inside the unit circle. 

A different method is based on the zeros of Pjq{y). When we have computed the 
N — 1 zeros y n of Pn(v), we can compute the corresponding zeros z n of Qn(z), by using 
the relation z n + l/z n = 2 — 4y n . We have (recall that Qjv(1) = 1) 

Pn(v) = Qn(z)Q n (1/z) = J] X - Y ^j Z 1 " ZZn 



-a 1 ~n 

n = l 

j-j 1 z n (z n + l/z n - z - 1/z) 



(1 z n f 

JV-1 

<*■>"- n ^ 



For each value y„ there are two z— zeros: z n and l/2„. We use the zeros z n inside the 
unit circle. By expanding the product representation of Qn we obtain the coefficients 
f(n) of (1.12), which are needed in (1.9). By expanding and using (1.7) and (1.9) we 
obtain the coefficients h(n) from 



1 /1 _l_ 1 / v v 



v n = v 7 ra = 

We verify this second method by taking N = 2. We have P2G/) = 1 + 2y, with zero 
?/i = — The corresponding z\ inside the unit circle is z\ = 2 — \/3- Equation (1.17) 
reads 

This gives the coefficients /i(n) given in (1.16). 

In Figure 1.1 we show the 99 zeros of Pwo(y) and the corresponding z— zeros. 

In Shen & Strang (1995) the location of the zeros of Pn(u) for large values of 
N is discussed. The zeros approach a limiting curve |4y(l — y)\ = 1 in the complex 
plane. As remarked by Shen and Strang, the wide dynamic range in the coefficients of 
Pn(u) makes the zeros difficult to compute for large N. They give pictures of the zeros 
located near the limiting curve up to N = 70, and they observe false results when using 
standard Matlab methods for N = 100. In the present paper we give an asymptotic 
expansion of the zeros, and we show by numerical verification that these approximations 
are excellent starting values for obtaining iterated high precision values of the zeros. The 
methods of computing the high degree polynomial Pn(v) are also based on a uniform 
asymptotic representation of the incomplete beta function, of which details are given in 
a final section. 
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Figure 1.1. Zeros y n (black dots) of Pn(v) for N = 100 and correspond- 
ing zeros z n (open dots) denned by z n + l/z n = 2 — 4y„. 

2. Asymptotic inversion of the incomplete beta function 

We recall the relation given in (1.15) 

P N (y) = (l-y)- p I 1 . y ( P ,p), (2.1) 

where the incomplete beta function is defined in (1.14). In Temme (1992) the asymp- 
totic inversion of the incomplete beta function I x (a,b) for large values of a and b is 
discussed. That is, an asymptotic approximation is given of the solution x of the equa- 
tion I x (a,b) = q, with < q < 1. We can use the same inversion method for finding 
complex a;— zeros (different from zero) by taking q = 0. In Temme (1992) three cases 
have been considered: 

• a and b large, with a — b bounded; 

• a and b large, with a/b and b/a bounded away from zero; 

• at least one of the parameters a, b is large. 

The first case applies to the present situation, because we are interested in the zeros of 

h-y(N,N). 

We repeat the main steps of the inversion procedure of I x (a, a + (3) given in §2 of 
Temme (1992), by taking a = N, (3 = 0. First we write h-yiN, N) in terms of an error 
function. 
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2.1. Transformation into Gaussian form 

We have 



'1-2/ 



a-N rl-y j f 



B(N,N)J L v n t(i-ty 
We transform this to a standard form with a Gaussian character by writing 

-|C 2 = ln[4i(l-i)], < * < 1, sign(C) = sign(i - ±), 
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In [4^/(1 - y)], < 2/ < 1, sign(7?) = sign(± - 2/). 



(2.2) 



Therefore 



-JV 



h-y(N,N) 



B(N, N) 

We can invert the relations in (2.2): 



y e " 



t(l - t) d( 



*=J l±v/l-exp(-iC 2 ) 



l± A /l-exp(-ir ? 2 ) 



l + C^/[l-exp(-iC 2 )]/C 2 _ 
1 - ??^/[l - exp(-i?7 2 )]/?7 2 



(2.3) 



where the second square roots are non-negative for real values of their arguments. It 
easily follows that 

1 dt_ = -C 
t(l -t) d( 1 - 2t' 

and that the following standard form (in the sense of Temme (1982)) can be obtained 



where 



We have 



$(iV) 



$(iV) 



i r(iv + i; 



Viv r(jv) 



HO 



1/2 
2^ 



- 1 "i iV " 1 + Tii iV " 2 + --- 



l-exp(-K 2 ) 



(N -» oo). 



(2.4) 



(2.5) 



(2.6) 



The function 0(£) is analytic in a strip containing 1R; the singularities nearest to the 
origin occur at ±2y^r exp(±i7r/4). The first part of the Taylor expansion of this even 
function reads: 

We write (2.4) in the form (see Temme (1982)) 

h- v (N,N) = ierfc (-riy/N/2^ + R N (v), 
where r\ is defined in (2.2), and erfc is the error function defined by 

erfc z = —= [ • 
V 71 " Jz 



dt. 



(2.7) 



(2.8) 
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In fact, we replace in (2.4) $(iV)0(£) by 1, and the error is contained in Rn(t)). We try 
to find zeros of 7 1 _ 1( (iV, iV) by using representation (2.7), assuming that N is large. First 
we find zeros in terms of tj; afterwards, we determine y from the second line in (2.3). 
When N is large, we consider the zeros of the error function as a first approximation to 
the zeros of the incomplete beta function. 

2.2. Inversion of the incomplete beta function 

Let 770 , 77 solve the equations 

Ierfc(-7 ?0A /pv) =0, 

( 2 - 9 ) 

ferfcj-TyfiVj +R N ( V ) = 0. 

We assume that 770 is known and that, for large N, this value is a first approximation 
for the value 77 that satisfies the second equation in (2.9). The relation between 77 and 
770 is written as 

V = Vo + s, (2.10) 
and we try to determine e. We can expand the quantity e in the form 

as iV — > 00. The coefficients can be found by using a perturbation method. For 
details we refer to Temme (1992). 
The first coefficient e\ is given by 

ei = - 1110(770), (2.12) 
^0 

where <j> is defined in (2.5). The quantity £\ is an odd function of 770 and analytic in the 
disc I r/o I < 2v / tt- Further terms are given by 

£ 2 = -^r(80ei+80'ei - - 40e?); (2.13) 

e 3 = ^^(i 28 ^ + 128^'eiei - 16^ei + 128<£'e 2 + 640"e?+ ^ 
- 128^'ei + - 1280ei£ 2 - Mtpelrj 2 - 64<ps 2 v4 ~ l&<psf). 

The derivatives <j>',e', etc., are with respect to jj, and all quantities are evaluated at 770- 
More information on the coefficients E{ will be given in §4.2. 

3. More details on the zeros and numerical examples 

As explained in the previous section, the zeros of Pn(v) are approximated in terms of 
the zeros of the error function. From Fettis, Caslin & Cramer (1973) we know 
that two infinite strings of zeros of erfc w occur in the neighbourhood of the diagonals 
v = ±u in the left half plane u < 0, w = u + iv. The first few zeros are given in Table 
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k 


Uk 




-i 
1 


— 1.354ol. 


.. ±« i.yyi4o... 


2 


-2.17704. 


.. ±i 2.69114... 


3 


-2.78438. 


.. ±i 3.23533... 


4 


-3.28741. 


.. ±i 3.69730... 


5 


-3.72594. 


.. ±i 4.10610... 



Table 3.1. First five pairs w k = Uk ± ivk of zeros of erfc«? 

3.1. Numerical values of the first 100 zeros of erfc w and asymptotic approximations of 
the zeros are also given by Fettis et al. A first order approximation reads 

wf ~ ^2ir(k - 1/8) e ±37 ™ /4 , k -» oo. (3.1) 

When ?7 is a zero of the right-hand side of (2.7), we use (2.10)— (2.11) for obtaining 
an asymptotic expansion of t] for large positive values of N; —t]o^/N/2 is a zero of the 
error function. The mapping of the t]— plane to the y— plane is given by (see (2.2) - 
(2.3)) 

-\r, 2 = In [42/(1 -y)]. (3.2) 

Because the zeros in terms of r\ occur in the neighbourhood of the diagonals 3ft?? = ±9?j 
in the right half of the r;— plane, we see that the y— zeros of Pn(v) occur near the curve 
defined by |4y(l — y)\ = 1, with Ufa/ < |. The full curve is a lemniscate; the extreme 
points cut the real axis at y = ^ ± 

The mapping in (3.2) is singular at the points 

±2Vk^e ± ^ n \ A; = 1,2,3,... 

on the diagonals in the rj— plane. All these branch points are mapped to y = Of 
special interest are the points on the diagonals in the right half plane, for instance the 
points 2y / 7r exp(±7r«/4). The parts of the diagonals given by n = pe ±vi / 4 ,0 <p< 20F 
are mapped to the left leave of the lemniscate, the parts satisfying 2y / 7r < p < 2v / 2~7r , 
are mapped to the right leave, and so on. 

If N is even Pjv(y) has one real zero (recall that this polynomial is of degree 
N — 1), and this real zero will be approximated by using w^^ 2 (see (3.1)), and for the 
complex zeros we can use w^, k = 1,2,... A^/2 — 1. When N is odd we use wf, k = 
l,2,...(JV-l)/2. 

We have used the asymptotic expansion (2.11) for N = 100 with five terms, and 
obtained approximations of the zeros of Pn(v) which were accurate with a precision of 
at least 12 digits. The zeros near y = \ are more accurate than those near the extreme 
point V = \ — The latter correspond to rj— values that are closer to the singular 
points 2y / 7r exp(±7r«/4) of the mapping given in (2.2). The asymptotic expansion (2.11) 
breaks down at these points. 

When N increases the used rj— zeros are bounded away from these singular points. 
To see this, observe that we need w^, with 1 < k < N/2. From (3.1) we see that w^, 2 ~ 
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VNtt e ±3 W 4 . Because wj., a zero of the error function, corresponds to —rjo^N/2 , we 
see that the maximal i]o satisfies r)o ~ y/2~Tr exp(±7r«/4). 



Remark 3.1. We will not discuss the role of other zeros in the rj— plane in the neigh- 
borhood of the diagonals. A similar phenomenon is explained in Temme (1995). The 
point is that, although we started with integer values of N, the incomplete beta func- 
tion in (2.1) and (2.7) can be interpreted with general complex values of N. When N 
is not an integer there are more y— zeros, but those occur on other Riemann sheets of 
the multi- valued incomplete beta function. 

3.1. Numerical aspects 

When N is large the representation of Pn(v) hi (1-1) is n °t suitable for verifying the 
accuracy of the zeros that are obtained by asymptotic methods with, say, 12 accurate 
digits. When N = 100 the coefficient of the power y" is an integer with about 60 
digits, and ( 1 gg)y 99 is about 10 30 when y = \. It follows that the zeros have to be very 
accurate to verify them by using (1.1), and we have to perform the computations with 
many extra digits. 

The incomplete beta function can computed for large values of N by using a uniform 
asymptotic expansion, which a gives a more stable representation. More details are given 
in the next section. We have verified the accuracy of the zeros by using this expansion, 
and obtained the accuracy mentioned above. However, we can use different numerical 
verifications, for instance the ones based on sums and products of zeros of polynomials. 
For the n zeros y k of the polynomial P(y) = a n y n + . . . + ao, we have the relations 



^ Uk = -<in-i, a n Y[ Dk = (-l) n «o- 



k=l k=l 

In the present case we know that the zeros y k of Pn(v) have to satisfy the rules 

1 f 2N ~ 2 \ TT\. r ^N-i 



^Vk 2 -, 

k = l " v " " 7 k = l 



U_Jn^=(-i) 

v 7 k=i 

We verified both rules for N = 100 and we obtained by using 5 terms in (2.11) 

1 -i o /198\ TT 

Vk + 2 = "2-56 • • • x lO" 13 , ( 9Q ) J] y k + 1 = -1.56 

k=l ^ 7 k=l 



x 10 -12 , 



which confirms our earlier claim about the accuracy of the zeros. 

We have obtained more accurate values of the zeros of P^(y), N = 100 by using the 
multi-length facilities of Maple and using the earlier obtained zeros as starting values in 
a Newton-Raphson process based on representations of Pn(v) given in (2.1) and (2.7). 
In this way we obtained the zeros y n with 45 relevant digits. 

In Shen &: Strang (1995) there is no discussion on numerical methods for ob- 
taining the filter coefficients h(n), once the zeros y n or z n are available. As we have 
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explained in §1.3 we need the f(n) in (1.17), but the computation of f(n) by using the 
zeros of Qn( z ) may be unstable when N is large. Another source of errors comes from 
expanding the product of the binomial term and Qn(z) in (1.17). Unfortunately, we 
don't have available asymptotic methods for obtaining the filter coefficients (we have 
powerful methods for Pjy, but not for Qn)- Further numerical research is needed to 
compute h(n). 

In fact, highly accurate values of y n are needed in order to obtain accurate values 
of h(n) by using straightforward methods. Expanding the right-hand side of (1.17) in 
powers of 1/z, and using 45-digits accurate values of z n , we have obtained a set h(n) 
that satisfies 

2JV-1 

^ /i 2 (ra) - 1 = 5.42 x 10 -19 , N = 100, (3.3) 

n=0 



which shows that the dominant h(n) are accurate with about 20 digits accuracy. A 
test based on the vanishing moments given in (1.10) indicated that only for rather low 
k— values a precision as in (3.3) could be obtained. 

In Figure 3.1 we give an idea of the values of the filter coefficients; we give the 
values of — 10 log \h(n)\,n = 0, 1, . . . , 199. The largest value of the coefficients is h(20) = 
0.39910 




Figure 3.1. Filter coefficients h(n) of the Daubechies wavelet -D 2 oo; 
shown are the values of - 10 log \h(n)\,n = 0,1,..., 199. 
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4. Asymptotic expansions 
4.1. The incomplete beta function 

We summarize from Temme (1982), (1996) the details of an expansion for the incom- 
plete beta function that is based on formula (2.7). We can expand 



This expansion holds uniformly in a strip |Ss7y| < V27T — 8, where 8 is a small positive 
number. The width of the strip is determined by the singularities of the mapping given 
in (2.2). As mentioned earlier, the singularities are defined by the t]— values for which 
exp(-i? ? 2 ) = 1. The r)— zeros needed for Pn(u) all he within this strip. 
By differentiating (2.4) and (2.7) with respect to 77, substituting 

00 

k = 

(a general form of (2.6)) and (4.1), we obtain a set of recursion relations for the coeffi- 
cients B k : 

1 - 4>{r)) „ dasetofB k 

B = , rjB k+1 = c k+1 <p(r)), A: = 0,1,2,... 

rj M) 

with <j) given in (2.5). The next few B k are 

1 



B 1 = -£^(3^>r? 2 -8</> 3 + 8) 
1 



8?? 



B 2 = -—^(25(f)7] 4 - 240r/ 2 </> 3 + 384</> 5 - 384). 
In these formulas we have removed the derivatives of <f) by using 



d<t>{rj) = 4>(y) 
drj rj 



(4.2) 



which easily follows from (2.5). 

For verifying the numerical computations of the zeros we have used expansion (4.1) 
with terms up to and including B w /N 10 . 

4.2. The coefficients e k of (2.11) 

In (2.14) and (2.15) we have given the coefficients e 2 and e 3 . It is convenient to remove 
derivatives by using (4.2). In this way we obtain 

e 2 = 1(8 + 3r? 2 - 8f + 4r/ 3 £l - 8fe lV - 4 e 2 r/ 2 )/r? 3 , 



13 

e 3 = —(-40 + 5?? 4 + 56(f) 4 - 16(f) 2 + 8rj 2 - 38i] 2 (f) 2 - 16r?£i + l6^eW - 32(f) 2 e^ 3 
16 

+ 480 4 eir/ - 6?? 3 ei + 16</> 2 £i?? + 4? ? 5 e 1 - 8?? 4 £ 2 </> 2 + 16<£ 2 eV + 8e??? 3 - 8e^ 4 )/?? 5 . 
To avoid cancellations for small values of \tj\ we can use the expansions 

1 1 3 1 y 1 11 1 15 1 17 

£l ~ 8 V ~192 V + 92160^ " 23224320^ + 4954521600^ ~ 980995276800 V +mmmm 
1 5 3 7 5 1 - 407 9 

£2 = V 71 H 71 H 71 71 + . . . , 

128 ' 1536 ' 40960 ' 81920 ' 371589120 ' 
_ -5_ 11 3 63 5 _ 823 7 

£3 "1024^ ~ 24576^ + 327680 V ~ 165150720^ + '"' 

- ~ 21 37 3 179 5 

64 "32768^ + 65536 V + 5242880 V + '"' 
399 219 o 

£5 = 71 H 71 + . . . . 

262144 ' 2097152 ' 
All these expansions have radius of convergence 2^/^^ . 
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